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Abstract 

A new algorithm has been developed to compute low Mach Numbers supercritical fluid flows. The algorithm is applied 
using a finite volume method based on the SIMPLER algorithm. Its main advantages are to decrease significantly 
the CPU time, and the possibility for supercritical fluid flow modelisation to use other discretisation methods (such 
as spectral methods and/or finite differences) and other algorithms such as PISO or projection. It makes it possible 
to solve 3D problems within reasonable CPU time even when considering complex equations of state. The algorithm 
is given after first a brief description of the previously existing algorithm to solve for supercritical fluids. The side 
and bottom heated near critical carbon dioxide filled cavity problems are respectively solved and compared to the 
previously obtained results. 



1. Introduction 

In the last decade, numerous numerical works 
have been devoted to the modeling of supercritical 
fluids (SCF) |1)2|3|4)5)6|7] . especially near their gas- 
liquid critical point. In fact, such fluids exhibit quite 
unusual properties, behaving as highly expandable 
gases with liquid-like density. An important point 
which has been pointed out is the apparition of 
the so-called piston effect in a closed supercritical 
fluid cell heated on a wall, where this piston ef- 
fect acts like a fourth mode of transport of energy. 
As a matter of fact, Onuki and al [I] pointed out 
the thermodynamic importance of the adiabatic 
heating while a more detailed hydrodynamic mech- 
anism of thermalization was proposed by Zappoli 
and al [2]. Close to the sample wall, heat diffusion 
makes a thin hot boundary layer expand and com- 
press adiabatically the rest of the fluid. A spatially 
uniform heating of the bulk fluid occurs pQwith a 
thermalization which should process at the velocity 
of sound [2]. Simultaneously, some dedicated ex- 
periments |8|9|10| have evidenced the existence of 
minute (a few iim/s, see [10]) but really efficient 



flows at the border of the expanding diffusive layer 
and compressed bulk fluid, independently of the 
closed geometries under consideration [10J. Such 
an effect is accounted for into the transport equa- 
tions through source terms, which are non-linearly 
related between density, temperature and pressure, 
and a real equation of state of supercritical fluid. 

The various research who have been conducted 
while treating critical fluids have been done at 
very low mach numbers M a (where M a = flow ve- 
locity/speed of sound). Depending on the Mach 
number and on the transient or not character of the 
flow, depending on the type of fluid, ideal or very 
expandable as supercritical fluids are, the relative 
strength of the coupling of mass, momentum, en- 
ergy and equation of state can be very disparate 
and one cannot use an universal efficient method for 
all the range of Mach number. As a matter of facts 
the dynamic coupling with pressure as introduced 
in high velocity flows "propagates" by the equation 
of state to the work of pressure forces in the energy 
balance because of the strong pressure gradient. In 
the case of very low velocity flows of ideal gazes 
this term is vanishing and the energy balance is 
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mainly governed by heat diffusion while the mo- 
mentum balance is that of a non compressible fluid 
in the Boussinesq approximation. However when 
it came necessary to compute unsteady flows of 
very expandable fluids like near critical fluids, new 
problems arose. These fluids are characterized by 
diverging isothermal expansion coefficient, isother- 
mal compressibility, thermal conductivity and heat 
capacity at constant pressure, their exponent being 
such that the heat diffusion goes to zero when near- 
ing the critical point of the phase diagram. Their 
isentropic compressibility, linked to the velocity of 
sound which goes to zero very slowly, remains com- 
parable to that of ideal fluids. This is why the term 
"hyper compressible fluids" as often encountered for 
supercritical fluid can be misleading and it is better 
to say "hyper expandable fluids". This mean first 
that the coupling between velocity and pressure is 
still very weak in the momentum balance equation. 
However due to the diverging isothermal expansion 
coefficient, the equation of state is no longer a pas- 
sive link between density, pressure and temperature 
since it makes the density variations in heat diffu- 
sion layer to be several order of magnitude larger 
than the temperature ones. Accordingly, the term 
of the energy balance which represents the work of 
the pressure forces due to the deformation of a fluid 
element during its motion becomes prominent, as it 
is in high velocity ideal fluid flows. The consequence 
of this fact was the numerical difficulties encoun- 
tered during the first attempts to calculate heat 
propagation by thermoacoustic coupling in super- 
critical fluids. In supercritical fluids where we work 
in a closed geometries, the heat input generates a 
piston type effect in which the isentropic compres- 
sion due to the slow motion generated by this piston 
effect is modified by acoustic waves. These waves 
propagate back and forth many times, because they 
are reflected at the walls of the domain. In that 
problem, we have one length scale, let say domain 
size, and two times scales: the long time it takes the 
slow diffusive flow to travel one length scale and the 
short time it takes an acoustic waves to travel one 
length scale. Therefore, the numerical approaches 
have opted for two solutions to completely treat 
heat transfer phenomena: the first one considers 
the simulation in time of the order of acoustic time 
and remains appropriate to account for initial heat- 
ing period, while the second one accounts for time 
higher than the piston effect time up to diffusive 
time. In the first case, the full transport equations 
are taken without any approximation, whereas in 



the second case in order to have in the numerical 
algorithm time steps which are not drastically small 
(in the order of acoustic time), one has to resort to 
the low mach number filtering approximation [3]. 
The low mach number filtering decouples the den- 
sity from the dynamic pressure which is the pressure 
part dependent both on space and time |1 1) 12) 13] . 
The other pressure part refers as a thermodynamic 
pressure in the following text and it is a part de- 
pendent only in time, except when its accounts for 
a possible local contribution due to the hydrostatic 
pressure which plays a significant role in supercrit- 
ical fluids under gravity field, even for very small 
heights of the cell [5] . As a consequence the numer- 
ical algorithms chosen have been limited to finite 
volume methods which require excessive amount 
of CPU time due to the necessary iterative scheme 
for coupling all equations, one has to keep in mind 
that due to the singularities arising near the critical 
point direct methods of resolution have failed [7] . 

In this paper we present an algorithm which de- 
couples at each time step the energy equation and 
the equation of state from the momentum and mass 
conservation equations. This decoupling brings a 
substantial reduction of the CPU time enabling a 
more straightforward extension to 3D modeling. It 
also permits to use other methods than finite vol- 
umes, i.e. finite differences and spectral methods, 
enables one to use besides SIMPLE and SIMPLER 
algorithms[14j, more direct algorithms such as PISO 
[15] or a modified projection technique [16], and fi- 
nally gives the possibility to extend simulation us- 
ing any complex equation of state (EOS). The algo- 
rithm is also applicable to several situations encoun- 
tered in low mach numbers flows as for example in 
subsonic combustion problems, or, more generally, 
in flows with velocities much smaller than speed of 
sound but having important density variation with 
temperature. 

The paper is organized as follows: Section 2 de- 
scribes the previous algorithm used in the model- 
ing of SCF under the low Mach number approxima- 
tion and then describe the new algorithm. Section 
3 shows numerical validation for two typical exam- 
ples, before concluding in Section 4. 

2. Numerical algorithm and Low mach 
number filtering 

In the present case, where we interest ourselves to 
critical fluids, we have to point out first the major 
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difficulties which arise when approaching the critical 
point. Going to the limit close to the critical point, 
the main physical properties are diverging, terms 
such as the pressure work in the energy equation 
are becoming as in compressible flows most leading 
terms and are in fact the terms in the equations re- 
sponsible for the "piston" effect with a uniform rise 
of the temperature in the bulk for closed domainspQ. 
In addition, the most straightforward equation of 
state which we can use to describe these real flu- 
ids is the van der Waals equation of state which is 
non linear and for which a direct resolution is not 
possible close to the critical point due to the fact 
that this equation degenerates close to the critical 
point [7]. Therefore one has to use iterative method 
to solve for the density, the thermodynamic pres- 
sure and the temperature. The convergence for this 
triplet (P — p — T) becomes quite slow due to the 
above and due to the fact also that in order to ob- 
tain the pressure work term in the energy equation, 
we have to obtain the divergence of the velocity field 
through the solving of the momentum equation and 
the equation of conservation of mass. 

Most of the interesting physics contained in study- 
ing the SCF are inherently transient, this brings us 
to another aspect which is: even tough we are in an 
hyper expandable fluid, the fluid velocities are much 
smaller than the speed of sound for the cases we 
would like to consider, i.e. small temperature pertur- 
bations of such fluids. Therefore another difficulty 
appears which this time puts the stringed on the size 
of the time step, the fluids being considered as low 
Mach number fluids [TT|12|13j . We then encounter 
the same difficulties as in low Mach number com- 
bustion, where the acoustic pressure waves force the 
algorithm for not becoming unstable to adopt small 
time steps: these time steps have to satisfy the CFL 
condition that requires a time step size smaller than 
the grid size times the reciprocal of the largest wave 
speed: 

, , 

At < - — 1 

max(c + \v\) 

where c is the speed of sound, v the flow velocity 
and Ax the grid size. 

There are two main approaches for solving low 
Mach number flows; the first one, is to use com- 
pressible solvers (density based) [17]; and the second 
one, is in extending incompressible solvers (pressure 
based) towards this regime[18|. Both of these tech- 
niques will anyway suffer from the pressure acous- 
tic waves, in order to alleviate these restrictions on 
time step two distinct techniques have been pro- 



posed, preconditioning and asvmptotic |17119)20|21| . 
Preconditioning techniques have many drawbacks 
for transient low Mach number flows and therefore 
have not been considered in this work. In the asymp- 
totic technique or perturbation approach[12j, a fil- 
tered form of the equations is employed to eliminate 
system stiffness. We expand all the variables in Tay- 
lor series in power terms of the Mach number, to do 
that we assume that the low Mach number asymp- 
totic analysis is a regular perturbation problem, i.e. 
all flow variables can be expanded in power series 
of M a (flow velocity / speed of sound) as for example 
the pressure which can be written as follow: 

P(X, t, M a ) = P Q (X, t) + M aPl (X, t) (2) 

+M 2 p 2 (X,t)+0(M!) 

Po, pi and p 2 are called the zeroth- (or leading), 
first- and second-order pressure respectively. 

The asymptotic analysis of the momentum equa- 
tion implies that 

VP = (3) 

Vpi = (4) 

and it shows that for the purpose of solving the fil- 
tered transport equations it is necessary only to re- 
tain the second order term in the above expansion 
([2]) which yield to: 

P(X, t, M a ) = P (t) + M 2 p 2 (X, t) + 0{Ml) (5) 

Equations ([3]) and |(4]) express that the zeroth-and 
first order pressure terms in the series expansion of 
the total pressure are independent on space and only 
dependent on time. 

The other variables which are T the temperature, 
p the density and V the vector velocity with compo- 
nents Ui, are expanded in the same manner in terms 
of M a number: 

V(X,t,M a ) = V (X,t) +M a V 1 (X,t) 

+M 2 V 2 (X,t) (6) 
+0(M 3 a ) 

T(X, t, M a ) = T (X, t) +M a T 1 {X, t) 

+M 2 a T 2 {X,t) (7) 
+0{Ml) 

p(X,t,M a )=p (X,t) +Mp al (X,t) 

+M 2 p 2 (X,t) (8) 
+0{Ml) 
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Again the asymptotic analysis shows that only the 
zeroth order terms have to be retained in equations 

Hi. 

Using the above asymptotic expansions and re- 
placing them in the transport equations of real flu- 
ids with van der Waals equation of state, we obtain 
the following dimensional equations (where g is the 
gravity vector ) : 

i) Mass conservation equation: 



dp 
dt 



dpU, 



dxi 

ii) Momentum equation: 



= 



(9) 



dpUi dpUj Ui 



Of 



djp2 
dxi 




( 9Ui 
\dx 3 



+ 



dU 3 
dxi 



(10) 



dxi 



pa 



iii) Energy equation (written here in term of C v 
and a viscous dissipation term <f>, see below): 



dpC v T dpUjC v T 
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dt 



dxi 



dxi 



in) 



iv) The van der Waals equation of state (with a 
and b as fluid-dependent parameters) : 



ap 2 = 



I -bp 



(12) 



The van der Waals equation leads to consider the 
isochoric specific heat Cy to be constant. Jointly 
to the use of this above classical equation of state, 
the singular behavior of the thermal conductivity A 
(estimated at constant critical density) is approxi- 
mated under the mean-field theory as follows: 

A (T) = Af, (T) + x MF T ^ 1 / 2 (13) 

where the label MF recalls for the mean-field value 
\ for the critical exponent of the power law in terms 
of t = T j Tc [T (T c ) is the (critical) temperature]. Ab 
is a background estimated far away from the critical 
point. The viscosity coefficient p (estimated at con- 
stant critical density) is approximated by its back- 
ground contribution in the mean field approxima- 
tion. Therefore, the viscous dissipation is written as: 



2// 



dU, 
dxi 

m, 

dx, 



dUj 
dx-j 



dxi 



dU, 
dxi 



More generally, we note that the van der Waals 
equation for thermodynamic properties as mean 
field approximations for transport properties does 
not only give us a phenomenological singular behav- 
ior for unusual properties such as compressibility, 
heat capacity at constant pressure, etc. but also 
needs to introduce analytical singular expression for 
the transport properties as for example for the ther- 
mal conductivity lfl3|) . However, we have considered 
the isochoric specific heat at constant volume Cy 
to be constant, as well as the viscosity p, because, 
for both of these two properties, their values are 
deduced from background contribution due to the 
fact that their critical exponents of divergence are 
very low and that they can be neglected for the 
values of distances to the critical point we are con- 
sidering in our numerical work (r > 1CP 4 ). Even 
tough the speed of sound is decreasing when com- 
ing close to the critical point, it stays however finite 
with values not lower than 50 ms -1 , implying that 
the Mach number stays of the order of 10~ 5 . There- 
fore under the approximation of low Mach number 
where the filtering of acoustic waves means that 
the time step in the numerical algorithm is bound 
only by the flow velocity as opposed to the speed 
of sound, the above supercritical fluid equations 
have been modeled exclusively with finite volume 
iterative methods [14]. The choice for such iterative 
methods was quite natural due to the fact that the 
main physical properties have strong temperature 
variations, and that the divergence term in the en- 
ergy equation was the main term in exhibiting the 
so-called "piston" effect. In supercritical fluids, we 
cannot avoid the iterative process coupling the ther- 
modynamic pressure, the temperature and the den- 
sity. This has implied that other numerical meth- 
ods such as spectral methods|l6], finite differences 
and schemes such as PISO or projection methods 
have been disregarded. In this paper, we will start 
from the previously existing numerical algorithm 
[3 J described below to introduce a new algorithm 
which takes fully in account the low Mach number 
nature of the flow in closed cavities for supercritical 
fluids, and through an Adams-Bashforth scheme 
|16|24 25 26 27j enables us to avoid the iterative 
process for the momentum and continuity equation 
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leading to a substantial saving in CPU times. 



2.1. Description of algorithm 1 

We have first developed Cartesian and Polar 
three dimensional codes where the discretisation 
of equations 



f9]fT2|) used a finite volume technique 
[14]leading to a finite set of algebraic equations that 
can be solved iteratively in a segregated manner. 
The solver chosen to accomplish such task are the 
BICGSTAB for temperature and velocities, and a 
preconditioned conjugate gradient for pressure. To 
couple momentum equation to continuity equation 
several methods are available, but when turning to 
finite volume, the most common methods used are 
the SIMPLE family type and its derivatives (SIM- 
PLEST, SIMPLEC, SIMPLER, PISO). The tem- 
poral scheme for convective part has several choices 
which are: hybrid scheme, power law scheme, Quick 
family schemes. For the examples shown later in 
this paper, they have been treated using the power 
law scheme. The SCF equations are inherently tran- 
sient and imply to choose an algorithm which will 
optimize the resolution at each time step. The ideal 
candidate should be a modified PISO algorithm (a 
predictor - corrector algorithm close in essence to 
the projection method of Chorin [21|16| ) where one 
does not need to iterate at each time step. Unfortu- 
nately as said earlier, the strong coupling between 
energy, density and thermodynamic pressure needs 
an iterative scheme and PISO cannot cope with 
the SCF system of equations for reasonable time 
steps without an outer iterative loop. In our case, 
the SIMPLE and SIMPLER algorithms have been 
chosen to treat the equations and their algorithm is 
rigorously the same as the original algorithms de- 
scribed by Patankar[14], the only special treatment 
being the solution of the density and the thermody- 
namic pressure which is shown hereafter. 

The van der Waals equation of state being non 
linear, if written as F(p) = 0, it has vanishing first 
derivative F'(p) near the critical point. The latter 
does not allow for a direct solution using methods 
such as Newton-Raphson. So we have to resort to 
a linearization of the equation of sate and solve the 
density iteratively. This done through the following 
linearization: 



= f(p k ) 



(14) 



or 



F( P ) 



F{ P ) 



(P + ap 2 )(l-bp) 



(Pa + ap 2 



(15) 



(16) 



T + b{P + ap 2 ) 

The convergence rate of such linearization can be 
found in Accary & Raspo [7] . 

In order to close the system, we need an equa- 
tion to derive the thermodynamic pressure Po- This 
equation is obtained through the conservation of 
mass as follow: 



pd£t 



Podfl 



(17) 



fl being the fluid domain and po the initial density. It 
has been found that internal iterations on equation 
(fl5j) or ([16| inside the SIMPLER iteration improve 
drastically the convergence stability and enables to 
take bigger time steps than in the case without inner 
iterations. When solving Po from equation lfl7|l and 
using Eq. l[T6|). Po in the numerator is at iteration 
k+1 and Po in the denominator at iteration k. 

We can summarize the different steps of the reso- 
lution as follow for each time step: 

(i) Solve the density field and the thermodynamic 
pressure using a known temperature 

(ii) Solve the momentum equations applying SIM- 
PLER or SIMPLE algorithm 

(iii) Solve conservation of mass 

(iv) Construct the divergence term and the viscous 
dissipation term going to the energy equation 
through the velocity just calculated at step 3 

(v) Solve the energy equation with the previously 
computed thermodynamic pressure at step 1 

(vi) Repeat to step 1 or achieve convergence on a 
convergence criteria calculated for each depen- 
dent variable 

The slowness of the convergence of the triplet 
(T, Po, p) impacts also on the momentum and 
mass conservation equations. It is to speed up con- 
sequently the convergence that we have adopted a 
new algorithm which is described here after. 

2.2. Description of the improved algorithm 2 

We will solve first the temperature equation (jTTJ) 
coupled with the van der Waals equation of state 
(fl"2|) and the divergence of the velocity. In order to 
obtain the divergence of the velocity, we will make 
use of the following: 



With 
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v.v = 



(1 - bp) 



dPg 
dt 



V. (AVT) - cj> 



- (P + ap 2 ) + 2ap 2 (1 - bp) (18) 

The above relation is obtained by taking the ma- 
terial derivative of the equation of state, and then 
substituting the appropriate terms from the mass 
and energy conservation equations, more details can 
be found later in the text equations (|21H25|) . Such a 
formulation of the divergence has been used to solve 
for low-frequency vibrations in a near critical fluid 
and transform the term source in the energy equa- 
tion without however changing the algorithm of res- 
olution |22|23j . 

We can notice that all the terms in Eq. (fl8|) can be 
calculated using only Po, p, and T at the exception 
of the viscous dissipation term In the low Mach 
number approximation the viscous term is of second 
order in terms of Mach number and can be neglected 
in our case, but for sake of generality, we will leave 
it in the equations. As described in the previous al- 
gorithm 1, the thermodynamic pressure Po can be 
obtained through the conservation of mass Eq. (fl7|) , 
the density through the van der Waals equation of 
state Eq. lfl2|) . If we want to have an efficient algo- 
rithm we have to decouple the momentum and mass 
conservation equations from the state and energy 
equation. In order to do that we will resort to an 
Adams-Bashforth scheme to linearize the (pV) n+1 
term in the - convective term and the dissipation 
term in the energy equation as follow: 



(pv) n+i = l(pv) n -U P vr 



(19) 

(20) 



2 yr ' 2 

V - -< 

T 2 s 

By linearizing only the (pV) n+1 instead of the full 
convective term including temperature we do not 
have to change the schemes previously used for the 
convective part (hybrid, power law, quick, smart 
etc.). The latter is specially useful for finite vol- 
ume methods. Whereas in the spectral methods we 
can use the Adams-Bashforth discretisation for the 
full advective-convective term, the stability of such 
a scheme is discussed in Ouazzani & al[24]. 

We have now decoupled at each time step the mo- 
mentum and mass conservation equations from the 
energy and state equations. 

It means that through step 1 to 6 at each time 
step the temperature can be solved independently of 
velocities at time n+1 but rather by using velocities 
at time n, n-1. 

The new algorithm will be as follow for each time 
step: 



Initialisation of variables 



Start of time step 



Construct the viscous dissipation term and the velocity 
for convective term using Adams-Bashforth 



Start of Temperature iteration 



Solve density and thermodynamic pressure 



Construct the divergence of velocity 



Solve energy equation 




Solve momentum and mass equations using 
SIMPLER 



Figure 1. Algorithm 2 organigram 

(i) Solve the density field and the thermodynamic 
pressure using a known temperature 

(ii) Construct the divergence term using equa- 
tion (fLS)) and the viscous dissipation term 
described above. 

(iii) Solve the energy equation with the previously 
computed thermodynamic pressure at step 1, 
and use equation (fl9|) for velocities in the con- 
vective term. 

(iv) Repeat to step 1 until convergence is achieved 
on temperature, thermodynamic pressure and 
density 

(v) if convergence is obtained go to step 6 

(vi) Solve the Navier-Stokes equations applying 
SIMPLER or SIMPLE algorithm (as shown 
in figure 1) 

The convergence speed of these equations be- 



6 



comes similar to the incompressible equivalent. 

By doing in such a manner, we have a better accu- 
racy in solving temperature, thermodynamic pres- 
sure and density; the piston effect can be shown 
without having to solve the momentum equation for 
time of the order of piston time. This approach en- 
ables us after optimization to solve more easily three 
dimensional problems. 

In a future paper, we will present a spectral 
method resolution of SCF and an extension to 3D 
cylindrical problems. 

This treatment can be applied similarly with other 
equations of state as follow: 

Let's assume that we have the following relation 
for the EOS: 

F(p,P,T) = (21) 
We can then derive this equation and obtain: 

dp OF dTdF dPdF 

— 1 1 =0 (22) 

dt dp dt dT dt dP y ' 

Then using the following two equations in equation 

(El): 

§--!(§£) V.V + Iv.(AVT, + * ,23) 
dt p \oT J p p p 



(24) 



We obtain an equation for the divergence of velocity: 



V.V 



gff + j£(v.(Avr)- 

O ( 9F , T (dP\ OF 
f \" dp + p \0T) p OT 



(25) 



It follows also if another EOS than van der Waals 
is chosen that the other physical transport proper- 
ties will have also to be redefined accordingly. 



in a square cavity of 1cm of side dimension filled 
with pure CO2 set at IK above the critical temper- 
ature. All its boundaries are thermally insulated ex- 
cept the one located at x=0 where the temperature 
is increased linearly of 10 mK over a period of Is. 
The results obtained with the two methods are very 
close (less than 0.1%) and compare also to the ones 
obtained by Zappoli & al [5]. 

The new algorithm has proved to be more effi- 
cient in terms of CPU time and a factor of 4 has 
been observed in the case of a mesh of 80 x 80. This 
speed up is also due to the fact that the resolution of 
the energy equation does not need to recompute at 
each iteration the "a„" coefficients appearing in the 
finite volume discretisation of the differential equa- 
tions. The iterative set coupling temperature, den- 
sity and pressure can even be improved by a Newton- 
Raphson type method. The optimization of this al- 
gorithm will be presented in a future work with a 
spectral code. 

In Fig 2, we present temperature field obtained 
at 4.5s. We will not discuss the physics related to 
the test cases due to the fact that they have already 
been discussed thoroughly in [3j. We can just point 
out that we observe a hot spot at the left corner 
of the cavity after the Is heating period and then 
it is convected for higher times which confirms the 
previous results obtained in [3] and by our work with 
algorithm 1 presented above. 

Through the examination of the steps in the algo- 
rithm, we can see that the piston effect is an adia- 
batic effect related to a thermodynamic relation be- 
tween the heat supplied to the system and the gener- 
ation of a pressure work resulting in a fluid motion. 
In the new algorithm, the effect is fully embodied 
in the energy equation coupled with an appropriate 
EOS and the conservation of mass at each time step. 



3. Numerical validation 

For numerical validation purpose, we have chosen 
two already existing cases which have been solved 
extensively by many authors. These two test cases 
are the adiabatic heated cavity from one side [5]and 
the 2D Rayleigh-Benard problem in a squared 
cavity [4|7l6] . 

3.1. The 2D adiabatic heated square cavity 



We consider the problem of the interaction be- 
tween gravitational convection and the piston effect 
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Figure 2. Isotemperature at time t=4.5s with algorithm 2. 
Distance to the critical point IK. 

3.2. The 2D Rayleigh Benard problem in a square 
cavity 

In this second test case, we consider the Rayleigh- 
Benard |4I7|6| problem which consists in heating 
from bottom a square cavity filled with a super- 
critical fluid. The fluid is maintained at a constant 
temperature in the upper boundary of the cavity 
whereas the side boundaries are adiabatic and im- 
permeable. 

The square cavity is a two dimensional cavity with 
10mm sides filled with CO2 on the critical isochore, 
initially at IK above the critical temperature. A 
mesh of 70 x 90 is considered with strong refinement 
of the mesh at the four walls. The time step is kept 
constant at 0.05. 

The results in this case agree well with those ob- 
tained by Amiroudine & al[4]. In figs 3 & 4, we show 
the temperature contours for an increase of lOmK 
of the bottom wall, for different times at 6.4s and 
8.5s, respectively. We observe the thermal plumes 
which appear at the hot and cold walls. The piston 
effect generates a cold boundary layer at the top of 
the cavity where gravitational instabilities develop 
as well as at the heated bottom wall. In figure 5, we 
can see the corresponding velocity field at time 8.5s. 

In this present case, we have the same speed up 
as for the side heated 2D cavity presented above 
but furthermore we can increase the time step by a 
factor of 10. which leads to a factor of 20 and more 



of CPU time gain. 
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>- 
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Figure 3. Isotemperature at time t=6.4s with algorithm 2. 
Distance to the critical point IK. Same results as with al- 
gorithm 1 with a timestep 10 times bigger. 




Figure 4. Isotemperature at time t=8.5s with algorithm 2. 
Distance to the critical point IK. 
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Figure 5. Velocity Field at time t=6s with algorithm 2. 
Distance to the critical point IK. 

4. Conclusion 

In this paper, we have developed a new algorithm 
to solve low Mach number flows. The algorithm 
takes advantage of the low Mach number filtering 
by computing the divergence of velocity from the 
three combined equations (equation of state, mass 
conservation and energy) , by doing so and using an 
Adams-Bashforth discretisation for the convective 
terms as well as for the viscous dissipation term, one 
can decouple at each time step the energy equation 
from the momentum and mass conservation equa- 
tions. This decoupling results in a better stability 
of the full algorithm and allows for a substantial 
saving of CPU time. Even if time step limitations 
are introduced through the CFL condition, these 
limitations do not lower the time step as compared 
to the fully implicit case which is in fact limited in 
the time step size by physical aspects. 

The other benefits of such an algorithm are to 
be able to use non iterative methods such as PISO, 
as well as other discretisation techniques such as 
pseudo spectral techniques, and it is easily applied 
to low mach number flows in general (as encountered 
in combustion problems) . Existing code for ideal gas 
law can be easily modified to tackle supercritical flu- 
ids and 3D problems can be tackled in a reasonable 
amount of CPU times. 

One other main reason to introduce such an algo- 
rithm for low Mach number critical fluids is to be 



able to treat problem with real equations of state 
other than van der Waals. These EOS are CPU time 
consuming and were not often used whereas now, we 
are introducing them in our 2D and 3D finite vol- 
ume code. 
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